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ABSTRACT 

Galaxies form in hierarchically assembling dark matter halos. With cosmological three dimensional 
adaptive mesh refinement simulations, we explore in detail the virialization of baryons in the concor- 
dance cosmology, including optically thin primordial gas cooling. We focus on early protogalaxies with 
virial temperatures of 10 4 K and their progenitors. Without cooling, virial heating occurs in shocks 
close to the virial radius for material falling in from voids. Material in dense filaments penetrates 
deeper to about half that radius. With cooling the virial shock position shrinks and also the filaments 
reach scales as small as a third the virial radius. The temperatures in protogalaxies found in adiabatic 
simulations decrease by a factor of two from the center and show fiat entropy cores. In cooling halos 
the gas reaches virial equilibrium with the dark matter potential through its turbulent velocities. We 
observe turbulent Mach numbers ranging from one to three in the cooling cases. This turbulence is 
driven by the large scale merging and interestingly remains supersonic in the centers of these early 
galaxies even in the absence of any feedback processes. The virial theorem is shown to approximately 
hold over 3 orders of magnitude in length scale with the turbulent pressure prevailing over the thermal 
energy. The turbulent velocity distributions are Maxwellian and by far dominate the small rotation 
velocities associated with the total angular momentum of the galaxies. Decomposing the velocity field 
using the Cauchy-Stokes theorem, we show that ample amounts of vorticity are present around shocks 
even at the very centers of these objects. In the cold flow regime of galaxy formation for halo masses 
below 1O 12 M0, this dominant role of virialization driven turbulence should play an important role in 
for star formation as well as the build up of early magnetic fields. 
Subject headings: Cosmology: high-redshift — galaxy formation — star formation 
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1. INTRODUCTION 

The process of virialization is clearly fundamental to all 
scales of galaxy formation. Lynden-Bell (1967) demon- 
strated that violent relaxation occurs during the virial- 
ization of a dissipationless system, but does the equiv- 
alent occur for the baryonic matter? If it does, how it 
achieves virial equilibrium should be inherently different 
because of hydrodynamical effects and radiative cooling. 
Additionally, this would create a Maxwellian velocity dis- 
tribution for the baryons as well. Turbulent velocities 
would exceed rotational ones. This would be at odds with 
the standard galaxy formation theories, which generally 
assume smooth solid body rotating gaseous distributions 
(e.g. Crampin & Hoyle 1964; Fall & Efstathiou 1980; Mo 
et al. 1998). The first occurrence of widespread star for- 
mation can be regarded as the commencement of galaxy 
formation, and its feedback on its host will affect all sub- 
sequent star formation. It is crucial to model the initial 
stage of galaxy formation accurately. Differences in ini- 
tial configurations of a collapsing halo may manifest itself 
in different types of central luminous objects, whether it 
be a stellar disk (Fall & Efstathiou 1980; Mo et al. 1998), 
a starburst (see §4 in Kennicutt 1998, for a review), or a 
massive black hole (Bromm & Loeb 2003; Volonteri et al. 
2005; Spaans & Silk 2006; Begelman et al. 2006). These 
differences may result from varying merger histories and 
the ensuing virial heating or turbulence generation of the 
new cosmological halo. 

For galaxy clusters, cosmological virialization has been 
studied extensively (Norman & Bryan 1999; Nagai & 
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Kravtsov 2003a; Nagai et al. 2003b; Schuecker et al. 2004; 
Dolag et al. 2005). It is customary to connect the veloc- 
ity dispersion to a temperature through the virial theo- 
rem for a collisionless system, where the potential energy 
equals twice the kinetic energy. However X-ray observa- 
tions and such cosmological simulations of galaxy clusters 
have indicated that turbulent energies are comparable to 
thermal energies. Central turbulent pressure decreases 
the density, but the temperature is largely unchanged. 
This leads to an increased entropy and a flatter entropy 
radial profile (Dolag et al. 2005) that is in better agree- 
ment with X-ray observations (e.g. Ponman et al. 1999). 
Simulations of merger dynamics suggest that turbulence 
is mostly generated in Kclvin-Helmholtz instabilities be- 
tween bulk flows and virialized gas during minor mergers 
(e.g. Ricker & Sarazin 2001; Takizawa 2005). Alterna- 
tively turbulence can be generated by conduction (Kim 
& Narayan 2003; Dolag et al. 2004) or acoustic transport 
of energy (Norman & Bryan 1999; Cen 2005). 

In standard galaxy formation models (Rees & Ostriker 
1977; Silk 1977; White & Rees 1978; Blumenthal et al. 
1984; White & Frenk 1991; Mo et al. 1998), gas shock- 
heats to the virial temperature as it falls into DM ha- 
los. These models succeed with considerable accuracy in 
matching various observables, such as star formation his- 
tories, galaxy luminosity functions, and the Tully-Fishcr 
relationship (White & Frenk 1991; Lacey & Silk 1991; 
Cole et al. 1994, 2000). Galaxy formation models depend 
on the virial temperature, most notably through the cool- 
ing function that controls star formation rates and their 
associated feedback mechanisms. Atomic hydrogen and 
helium radiative cooling is efficient in halos with masses 
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between 10 s and 10 12 M Q as cooling times can be less 
than the dynamical time of the system, a condition that 
galaxy clusters do not satisfy. This strong cooling sug- 
gests that in galaxies thermal energy may be less impor- 
tant for virialization than turbulent kinetic energy. Mo- 
tivated by the results of galaxy cluster turbulence, we 
investigate this potentially important role of turbulence 
and radiative cooling in galaxy formation, using a series 
of high resolution numerical simulations of protogalactic 
halos in this work. 

We consider kinetic energy and pressure forces in our 
virial analysis of protogalactic halos (cf. Shapiro et al. 
1999; Iliev & Shapiro 2001). This allows us to inves- 
tigate the equilibrium throughout the entire halo and 
determine the importance of each energy component in 
the virial theorem. Kinetic energy can be decomposed 
into radial and azimuthal motions along with turbulence, 
which affects the collapse of gas clouds primarily in three 
ways. First as seen in galactic molecular clouds, turbu- 
lence plays an integral part in current theories of star for- 
mation as the density enhancements provide a favorable 
environment for star formation (Larson 1979, 1981; My- 
ers 1999; Goldman 2000). Second if the turbulence is su- 
personic, gas dissipates kinetic energy through radiative 
cooling, which aids the gaseous collapse (Rees & Ostriker 
1977). Conversely turbulent pressure adds an additional 
force for the collapsing object to overcome and can delay 
the collapse into a luminous object. Last, turbulence pro- 
vides an excellent channel for angular momentum trans- 
port as the halo settles into rotational equilibrium to 
satisfy Rayleigh's inviscid rotational stability argument 
(Raylcigh 1920; Chandrasekhar 1961) in which the spe- 
cific angular momentum must increase with radius. Cos- 
mological hydrodynamic simulations have just begun to 
investigate angular momentum transport within turbu- 
lent collapsing objects, and turbulence seems to play a 
large role in segregating low (high) angular momentum 
gas to small (large) radii (Norman & Bryan 1999; Abel 
et al. 2002; Yoshida et al. 2006b). 

We study idealized cases of structure formation where 
stellar feedback is ignored because it provides a conve- 
nient problem to focus on the interplay between cos- 
mological merging, hydrodynamics, and cooling physics 
during the assembly of early halos. Some of the dis- 
cussed physical principles should, however, be applicable 
to galaxies of all masses. These simulations provide the 
simplest scenario to which we can incrementally consider 
further additional physics, such as H2 and HD cooling 
physics (Saslaw & Zipoy 1967; Palla et al. 1983; Flower 
et al. 2000), primordial stellar feedback (Whalen et al. 
2004; Kitayama et al. 2004; Alvarez et al. 2006; Yoshida 
et al. 2006a; Abel et al. 2007), metal enrichment from pri- 
mordial stars (Heger & Woosley 2002; Tumlinson 2006), 
AGN feedback (Springel et al. 2005; Kuhlen & Madau 
2005), and "normal" metal-enriched star formation (see 
Larson 2003, for a review). 

We present a suite of adaptive mesh refinement sim- 
ulations that are described in §2. Then we analyze the 
local virial equilibrium and shocks in halos in §3. There 
we also differentiate between infall through voids and fil- 
aments and its associated virialization. We discuss the 
situations in which virial heating and turbulence occur. 
Next in §4, we decompose the velocity distribution in 
principle axes to explore virialization in both the DM 



and baryonic components. Furthermore we decompose 
velocities into shear and compressible flows to study tur- 
bulent flows in the virialized gas. We discuss the impli- 
cations of these results on star and galaxy formation in 
§5. Finally we summarize in the last section. 

2. THE SIMULATIONS 

To investigate protogalactic (T v ; r > 10 4 K) halo virial- 
ization in the early universe, we utilize an Eulerian struc- 
tured, adaptive mesh refinement (AMR), cosmological 
hydrodynamical code, Enzo 1 (Bryan & Norman 1997, 
1999; O'Shea et al. 2004). Enzo solves the hydrodynam- 
ical equations using the second order accurate piecewise 
parabolic method (Woodward & Colella 1984; Bryan et 
al. 1995), while a Riemann solver ensures accurate shock 
capturing with minimal viscosity. Additionally Enzo uses 
an adaptive particle-mesh n-body method to calculate 
the dynamics of the collisionlcss dark matter particles 
(Couchman 1991). Regions of the simulation grid are 
refined by two when one or more of the following con- 
ditions are met: (1) Baryon density is greater than 3 
times flbpoN^ 1 ^, (2) DM density is greater than 3 

times ^cdmA)-^ 1+ ^\ and (3) the local Jeans length 
is less than 4 cell widths. Here N = 2 is the refinement 
factor; I is the AMR refinement level; <j> — —0.3 causes 
more frequent refinement with increasing AMR levels, i.e. 
super-Lagrangian behavior; p = 3Hq/8ttG is the critical 
density; and the Jeans length, Lj — y/l5kT/4irpGp,m H , 
where Ho, k, T, p, /1, and run are the Hubble con- 
stant, Boltzmann constant, temperature, gas density, 
mean molecular weight in units of the proton mass, and 
hydrogen mass, respectively. The Jeans length refine- 
ment ensures that we meet the Truelove criterion, which 
requires the Jeans length to be resolved by at least 4 cells 
on each axis (Truelove et al. 1997). 

We conduct the simulations within the concordance 
ACDM model with WMAP first year parameters 
(WMAP1) of h = 0.72, fl A = 0.73, M = 0.27, n b = 
0.024/i~ 2 , and a primordial scale invariant (n — 1) power 
spectrum with eg = 0.9 (Spergel et al. 2003). h is the 
Hubble parameter in units of 100 km s _1 Mpc -1 . Ct\, 
Qmj and Q& are the fractions of critical energy density of 
vacuum energy, total matter, and baryons, respectively. 
Last us is the rms of the density fluctuations inside a 
sphere of radius 8h~ x Mpc. 

Using the WMAP1 parameters versus the signifi- 
cantly different WMAP third year parameters (WMAP 3; 
Spergel et al. 2006) have no effect on the evolution of 
individual halos as are considered here. However these 
changes play an important role in statistical properties. 
For example, halos with mass 1O 6 M at redshift 20 cor- 
respond to 2.8(7 peaks with the WMAP1 but are 3.5u 
peaks for WMAP3. The Um/^ ratio also only changed 
from 6.03 to 5.70 in WMAP3. 

We also have verified that there is nothing atypical 
about the mass accretion rate histories of the objects 
we study. The mass accretion history of these objects 
exhibit smooth growth during minor mergers and accre- 
tion and dramatic increases when a major merger occurs. 
This behavior is consistent with typical halo assemblies in 
extended Press-Schechter calculations (Bond et al. 1991; 

1 See http://http://lca.ucsd.edu/portal/software/enzo 
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TABLE 1 
Simulation Parameters 



Name 


1 

[Mpc] 


z encZ 


N. 


part 


N grid 




N ceii 


Cooling model 


AO 


1.0 


15.87 


2.22 


X 


10 7 


30230 


9.31 


X 


10 7 (453 3 ) 


Adiabatic 


A6 


1.0 


15.87 


2.22 


X 


10 7 


40486 


1.20 


X 


10 8 (494 3 ) 


H,He 


A9 


1.0 


18.74 


2.22 


X 


10 7 


45919 


1.21 


X 


10 s (495 3 ) 


H,He,H 2 


BO 


1.5 


16.80 


1.26 


X 


10 7 


23227 


6.47 


X 


10 7 (402 3 ) 


Adiabatic 


B6 


1.5 


16.80 


1.26 


X 


10 7 


21409 


6.51 


X 


10 7 (402 3 ) 


H,He 


B9 


1.5 


23.07 


1.26 


X 


10 7 


20525 


5.59 


X 


10 7 (382 3 ) 


H,He,H 2 



Note. — Col. (1): Simulation name. Col. (2): Comoving box size. Col. (3): 
Final redshift. Col. (4): Number of dark matter particles. Col. (5): Number of AMR 
grids at the final redshift. Col. (6): Maximum number of unique AMR grid cells. 
Col. (7): Cooling model. 



Bower 1991; Lacey & Cole 1993, 1994; van den Bosch 
2002) and cosmological numerical simulations (e.g. Dc 
Lucia et al. 2004; Gao et al. 2005). The mass accretion 
histories in our simulations are well described by the fit- 
ting function of van den Bosch with Mq = 3 x 1O 7 M , 
Zf = 17, and v = 12.5. We also compare our data against 
the mass accretion histories of Gao et al., who tested 
their data against an extended Press-Schechter calcula- 
tion of the growth history of the halos. We find no major 
discrepancies between the two histories. 

The initial conditions of this simulation are well- 
established by the primordial temperature fluctuations in 
the cosmic microwave background (CMB) and big bang 
nucleosynthesis (BBN) (Hu & Dodelson 2002; Buries et 
al. 2001, and references therein). 

We perform two realizations with different box sizes 
and random phases. In the first simulation (simulation 
A), we set up a cosmological box with 1 comoving Mpc 
on a side, periodic boundary conditions, and a 128 3 top 
grid with three nested child grids of twice finer resolution 
each. The other simulation is similar but with a box 
side of 1.5 comoving Mpc (simulation B). We provide a 
summary of the simulation parameters in Table 1. These 
volumes are adequate to study halos of interest because 
the comoving number density of >10 K halos at z = 10 
is ^6 Mpc~ 3 according to an ellipsoidal variant of Prcss- 
Schcchter formalism (Sheth & Tormen 2002). We use the 
COSMICS package to calculate the initial conditions at z 
= 129 (119)* (Bertschinger 1995, 2001), which calculates 
the linear evolution of matter fluctuations. We first run 
a dark matter simulation to z = 10 and locate the DM 
halos using the HOP algorithm (Eisenstein & Hut 1998). 
We identify the first dark matter halo in the simulation 
with T v j r > 10 4 K and generate three levels of refined, 
nested initial conditions with a refinement factor of two, 
centered around the Lagrangian volume of the halo of 
interest. The nested grids that contain finer grids have 8 
cells between its boundary and its child grid. The finest 
grid has an equivalent resolution of a 1024 3 unigrid. This 
resolution results in a DM particle mass of 30 (101) M 
and an initial gas resolution of 6.2 (21) M@. 

Enzo employs a non-equilibrium chemistry model 
(Abel et al. 1997; Anninos et al. 1997). We conduct three 
simulations for each realization with (i) the adiabatic 
equation of state with an adiabatic index 7 = 5/3, (ii) 
a six species chemistry model (H, H + , He, He + , He ++ , 

t To simplify the discussion, simulation A will always be quoted 
first with the value from simulation B in parentheses. 



e ), and (iii) a nine species chemistry model that adds 
H2, H"|", and H~ to the six species model. In the nine 
species model, we use the molecular hydrogen cooling 
rates from Galli & Palla (1998). These models are differ- 
entiated in the text by denoting 0, 6, and 9, respectively, 
after the simulation name (e.g. simulation B0). Comp- 
ton cooling and heating of free electrons by the CMB 
and radiative losses from atomic and molecular cooling 
are also computed in the optically thin limit. 

To restrict the analysis to protogalactic halos in the 
H 2 models, we suppress H 2 formation in halos that can- 
not undergo Lya cooling by reducing the residual elec- 
tron fraction to 10~ 12 instead of a typical value of ~ 10~ 4 
only at the initial redshift (Shapiro et al. 1994). This 
mimics an extreme case where all H 2 is dissociated by an 
extremely large radiation background, and the halo can 
only collapse and form stars when free electrons from 
ionized hydrogen can catalyze H 2 formation. 

We end the simulations with non-equilibrium cooling 
when the gas begins to rapidly cool and collapse. We 
choose a final resolution limit of ^3000 (4000) proper 
AU, corresponding to a refinement level of 15. We end 
the adiabatic simulations at the same redshift. In a later 
paper, we will address the collapse of these halos to much 
smaller scales. 

3. VIRIAL ANALYSIS 

The equation of motion for an inviscid gas in tensor 
notation reads: 



Dvi 



d 



(1) 



where D / Dt = d/dt + Vjd/dxj is the total derivative. 
Here v is velocity; p is pressure; p is density; and g — 
V$ where $ is the gravitational potential. From this 
Chandrasekhar & Fermi (1953) derived the general virial 
theorem for a region contained within a surface S, in 
scalar form, 

1D 2 I 



2 Dt 2 



2T + V + 3(7-1)5- / pr-dS, (2) 



where 



V 



p( x )p( x '; 



dxdx' , 



(3) 



T = I J pv 2 dx, I = J px 2 dx, and £ = J e dx denote 
the gravitational potential energy, the trace of the ki- 
netic energy and inertia tensor, and the total internal 
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TABLE 2 
Halo Properties 



Name 


z coll 


M tot 


Pc 


T . a 
J vir 


T c 


<T> 






[Mq] 


[cm- 3 ] 


[K] 


[K] 


[K] 


H2 Induced Collapse 


AO 


18.74 


9.8 x 10 6 


8.1 


9200 


10000 


5700 


A6 


18.74 


9.8 x 10 6 


17 


9200 


7700 


5500 


A9 


18.74 


9.8 x 10 6 


1.6 x 10 6 


9200 


590 


5000 


BO 


23.07 


6.2 x 10 6 


13 


8300 


12000 


6000 


B6 


23.07 


6.2 x 10 6 


15 


8300 


8900 


5600 


B9 


23.07 


6.7 x 10 6 


3.0 x 10 6 


8700 


580 


4200 


Lye? Induced Collapse 


AO 


15.87 


3.6 x 10 7 


4.9 


19000 


17000 


12000 


A6 


15.87 


3.6 x 10 7 


1.8 x 10 6 


19000 


8700 


7300 


BO 


16.80 


3.5 x 10 7 


3.8 


19000 


31000 


11000 


B6 


16.80 


3.6 x 10 7 


4.0 x 10 6 


20000 


9000 


7500 



Note. — Col. (1): Simulation name. Col. (2): Rcdshift of 
collapse through H2 or Ly« cooling. Col. (3): Total mass. Col. (4): 
Central density. Col. (5): Virial temperature (i.e. eq. [6]). Col. 
(6): Central temperature. Col. (7): Mass-averaged temperature of 
the entire halo. 

a Virial temperatures arc calculated with = 1.22 in all cases. 

thermal energy, respectively. The surface term E s (the 
last term in eq. [2]) is often negligible in the outer re- 
gions of the halo. The system is not necessarily in virial 
equilibrium if I — 0, but the time-averaged quantity is 
zero when the entire system is in virial equilibrium. A 
system is expanding or contracting whether I is posi- 
tive or negative, respectively, based on energy arguments. 
Ballesteros-Paredes (2006) gives counterexamples to this 
simple interpretation. However, in the cases presented 
here, spherically averaged radial velocities are always 
negative. 

We define the halo as the material contained in a sphere 
with a radius r2oo enclosing an average DM overdensity 
of 200 and as such relates to mass by 

1/3 



^200 



GM 



(4) 



where M is the mass of the halo, Ocdm(z) is evaluated 
at a redshift z, and H is the Hubble parameter at z. The 
region where the cooling time is shorter than a Hubble 
time is denoted as the cooling radius r coo i, 

i C ooi(r C ooi) = H{zY x (5) 

(White & Frenk 1991). Mass and radius define a circular 
velocity and virial temperature, which are 



GM 
^200 



and T, r 



2k 



(6) 



for a singular isothermal sphere (see Bryan & Norman 
1998, with 13 = 1 and A c = 200). 0.59 and 1.22 in units 
of the proton mass arc appropriate values for p for the 
fully ionized and completely neutral states of a primor- 
dial hydrogen and helium mixture of gas, respectively. 
We use p — 1.22 throughout this paper. We note that 
Iliev & Shapiro (2001) considered non-singular, trun- 
cated isothermal spheres, and the resulting virial temper- 
ature is ~15% lower than the one calculated in equation 
(6). 

For 7 = 5/3, T v ; r is the temperature at which an ideal 
adiabatic gas reaches virial equilibrium with the specified 



potential. Please note that for an isothermal gas where 7 
is close to unity virial equilibrium is established between 
the turbulent energies and the gravitational potential as 
the 3(7 — l)£ term in equation (2) goes to zero. 

3.1. Local Analysis 

We evaluate the terms of equation (2) with respect to 
radius (i.e. the volume contained in a radius r). Figure 1 
illustrates the radial structure of (a) the turbulent Mach 
number, 



Mturb = 



pm h 



(7) 



(b) turbulent and (c) thermal energies per rrih, and (d) 
a "virialization" parameter 3 



a _ 3( 7 -l)£ + 2T 
' E s - V 



(8) 



of the adiabatic and radiative cooling simulations when 
the cooling halo collapses. Here v rm s is the three- 
dimensional rms velocity and is assessed using the gas 
velocities relative to the mean gas velocity of each spher- 
ical shell. In the top row of Figure 1, we also plot the 
Mach number, using v rms with respect to the mean ve- 
locity of the gas within r2oo- In the six and nine species 
simulation, this occurs at z c ,H2 = 18.7 (23.4) and z Ci Lya = 
15.9 (16.8), respectively. The radial profiles are centered 
on the densest point in the simulation with the collaps- 
ing halo. Several properties of the most massive halo in 
each simulation are detailed in Table 2. The sections of 
the Table compare the halo in the adiabatic, Lya, and 
H 2 simulations. 

Figure 2 shows the mass- weighted radial profiles of gas 
mass enclosed and gas density at z = z c ^ ya in the adia- 
batic and cooling cases. Both realizations are remarkably 
similar. Halos in the adiabatic case have a central core 
with a radius ~50 pc and gas density of ~5.0 (3.5) cm~ 3 . 
Core densities in simulation A are slightly higher than 
simulation B, which has larger thermal and turbulent 
pressures (see Figure 1). With radiative cooling, gas in- 
falls rapidly as it cools and undergoes a self-similar col- 
lapse with p oc r - 12 / 5 . 

3.1.1. Virial Radius 

We define the virial radius r vir when (3 = and 
dfi/dr < 0. When we radially average the Lya halo, 
r v i r = 419 (787) pc in the adiabatic simulations where 
the corresponding r2oo value is 615 (576) pc. 

A well defined shock exists between the interface be- 
tween voids and the halo. This material shock-heats to 
T v i r and virializes at a radius comparable to r 2 oo in the 
adiabatic cases. When we include radiative cooling, this 
radius decreases everywhere around the halo-void virial 
shock and is low as r 2 oo/2. In contrast to the voids, the 
filamentary gas shock-heats at an even smaller radius. 
Dekel & Birnboim (2006) also studied the stability of 
cold inflows within a hot virialized medium and found 
similar results. Figures 2 and 3 illustrate these changes. 

3 This is a modified version of the j3 used in Shaw et al. (2006) 
to account for kinetic energies and so that it does not diverge when 
V — > in the center. It still has the same behavior of f3 — > as 
j'->0. 
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Fig. 1. — A comparison of (top to bottom) turbulent Mach numbers, turbulent and thermal energies, and virial parameters between 
simulations with radiative cooling (dashed) and adiabatic models (solid). The main coolant is listed at the top of each column. The first 
and second columns display the state of these variables at z = z c ,H2 = 18.7 (23.1) for simulation A and B, respectively. The third and 
fourth columns are the data at z = z c L ya = 15.9 (16.8). The top row depicts the importance of radiative cooling in generating trans- and 
super-sonic turbulence throughout the halo during virialization. The thick lines represent the turbulent Mach number (eq. [7]), and the 
thin lines show the Mach number using the rms velocity with respect to the mean velocity of the halo. The middle two rows show that 
when radiative cooling is efficient the halo cannot virialize through heating but must virialize by increasing its kinetic (turbulent) energies. 
The dotted line in the third row marks T v i r (cq. [6]) with p = 1.22, estimated from the total halo mass. We plot the virialization parameter 
/3 (eq. [8] ) to investigate the local virial equilibrium (f3 = 0) , particularly at r2oo ■ Furthermore, j3 allows us to determine the mass-averaged 
dynamics of the system at a given radius, where j3 > and < correspond to decelerating and accelerating collapses, respectively. r200 (eq. 
[4]) is marked on the bottom of each column. 



At r200) densities in the adiabatic case begin to increase 
more rapidly than the cooling case as material accretes at 
the virial shock. No significant increase in dp/ dr is seen 
in the cooling case, indicative of a self-similar collapse. 

3.1.2. Adiabatic Model 

We start with the discussion of the adiabatic model as 
it is the simplest case and later compare the calculations 
with radiative cooling to this model. Virialization should 
transfer potential energy to kinetic energy that dissipates 
in shocks to thermal energy, which is the implication of 
the dissipationless virial theorem. The solid lines in Fig- 
ure 1 represent the energies in adiabatic models. The 
physics illustrated in this Figure are as follows: 

1 . Thermal energy — The gas shock- heats to the virial 
temperature at the virial shock. Virial heating contin- 
ues with decreasing radius as the surface term becomes 
significant in the interior. The resulting central tempera- 
ture of the halo is 10000 (12000) K, which is 1.2 (1.5) T vir , 
at a redshift of z Ci H2 when the H2 model collapses. At the 
time (z = z C! L yQ ) of collapse caused by Lya cooling, the 
central temperature is 17000 (31000) K, corresponding 
to 0.9 (1.6) T vir . 



2. Kinetic energy — It increases along with the thermal 
energy during virialization. The gas is generally turbu- 
lent, appearing as a velocity dispersion with a bulk radial 
inflow. At r2oo, the kinetic energy is equivalent to the 
thermal energy, T/£ ~ 1. This ratio steadily drops to- 
ward the center, where T/E ~ 1/3. This decrease in 
kinetic energy is apparent in all the calculations except 
simulation B at z C! Lya, increasing by a factor of two in 
the center. 

3. Turbulent Mach number — At T2oo, the turbulent 
Mach number Mturb lS maximal and varies from 1-3 in 
all adiabatic simulations. Mturb decreases to subsonic 
values ^0.15 but never below in the interior. Note that 
Mturb does not increase as the turbulent energy towards 
the center in simulation B because of the also growing 
sound speed there. 

4. Virialization parameter — Virial equilibrium is 
quantified by the virialization parameter /3, where the 
collapse is retarding or accelerating when it is negative 
or positive, respectively. At z Ci H2 and z C! Ly Q and in both 
simulations, (3 is within 20% of being virialized {(3 = 0). 
For comparison purposes, this corresponds to a halo hav- 
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Fig. 2. — Mass-weighted radial profiles for gas mass enclosed 
(top) and number density (bottom) for simulations A (solid black) 
and B (dashed red) at z = z Ci L ya = 15.9 (16.8). The light and heavy 
lines represent data for adiabatic and cooling models, respectively. 
The virial shock in the cooling halos occurs at ^r200, illustrated 
by the density increasing at smaller radii. We mark the radius T200 
= 615 (576) for simulation A (B). 



ing 80% of the required velocity dispersion for virializa- 
tion in the dissipationless case. At r v i r , (3 is nearly zero 
which defines the virialized object. Within r v - 1T , the val- 
ues decrease to values around -0.1 but stays < 0. 

Characteristics of turbulence in our adiabatic models 
are similar to ones found in galaxy cluster simulations 
(Norman & Bryan 1999; Dolag et al. 2005). Both groups 
find that turbulence provides ^5-30% of the total pres- 
sure, i.e. T/(T + £), in the cluster cores. Our proto- 
galactic halos have ~25% of the pressure in the turbu- 
lent form. Also the galaxy clusters in Norman & Bryan 
(1999) have comparable Mach numbers of ~1.6 at r v ; r , 
~0.5 at r v i r /3, and ~0.3 in the core. These similarities 
suggest that virial turbulence is generated over a large 
range of mass scales. 

3.1.3. Lya Cooling Model 

At z Ci Lya, halos in calculations with the H+He cool- 
ing model start to rapidly collapse. The dashed lines in 
the third and fourth columns of Figure 1 illustrate the 
energies of this model. 

1. Thermal energy — Compared to the adiabatic mod- 
els, the gas can radiatively cool through Lya emission 
to T ~ 8000 K within r coo i ~ r v - 1T . The entire halo is 
isothermal at this equilibrium temperature. Below this 
temperature, the cooling function of pristine gas drops by 
several orders of magnitude, and the gas can no longer 
cool efficiently. The thermal energy is ~65% lower than 
the adiabatic case. 



2. Kinetic energy — In response to the lesser thermal 
energy, the system tends toward virial equilibrium by 
increasing kinetic (turbulent) energy. The gravitational 
potential and surface terms do not appreciably change 
with the inclusion of radiative cooling. Turbulent en- 
ergy within r v j r increases as much as a factor of 5 when 
compared to the adiabatic case. 

3. Turbulent Mach number — The changes in thermal 
and kinetic energies equate to a increase of Adturb by a 
factor of 2-3 to values up to 1.5. The turbulence is super- 
sonic in all cases at the virial shock, but when we include 
radiative cooling, this trait emanates inward as the halo 
begins to rapidly cool. When the central core becomes 
gravitationally unstable, the entire halo is supersonically 
turbulent. 

4. Virialization parameter — The increased kinetic en- 
ergies compensate for the loss in thermal energy and the 
halo remains in a similar virial state. This is apparent in 
the remarkably similar radial characteristics of (3 in the 
adiabatic and H+He models of simulation A. 

3.1.4. H 2 Cooling Model 

The collapses caused by H2 cooling at z = z Ci H2 have 
very similar dynamics as the halos described in the pre- 
vious section. The dashed lines in the first and second 
columns of Figure 1 illustrate the energies of this model. 

1 . Thermal energy — H2 cooling is efficient down to 300 
K, so gas can depose a much larger fraction of its thermal 
energy. Inside r coo i ~ 0.32 (0.19) r v j r , thermal energies 
are only 5% of the values in the adiabatic models. 

2. Kinetic energy — The turbulent energies must in- 
crease as in the Lya case, and they increase by 93% 
(44%) on average inside r coo i. 

3. Turbulent Mach number — Similarly, Aiturb in- 
creases up to a factor of 10 to become supersonic at values 
up to 3 throughout the halo. They are somewhat larger 
than the Lya cases since H2 can cool to significantly 
lower temperatures than the virial temperature. 

4. Virialization parameter — The virial equilibrium of 
the halos are also similar to the other models. j3 smoothly 
transitions from nearly equilibrium at r v j r to an increased 
radial infall with f3 = -0.2 at 70 pc. Then it increases 
to 0.4 inside 10 pc, which corresponds to the gas decel- 
erating from the rapid infall as it encounters the central 
molecular cloud. 

3.1.5. Model Summary 

Baryons are close to virial equilibrium over three or- 
ders of magnitude in length scale by gaining both thermal 
and kinetic energies independent of cooling physics. Cen- 
tral temperatures of the adiabatic simulations are up to 
twice the nominal virial temperature. Similar to galaxy 
cluster studies, turbulence in the adiabatic model con- 
tributes ^25% to the energy budget with Mach numbers 
~ 0.3 in the center. In cooling cases, atomic and molecu- 
lar cooling inhibit virialization through heating, therefore 
the object must virialize by gaining kinetic energy up to 
five times the energy seen in the adiabatic models. This 
translates into the flow becoming supersonically turbu- 
lent with Mach numbers ranging from one to three. 
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Fig. 3. — Differences in entropy and density in a protogalactic halo at z = z c ,Lyce Left to right: simulations AO, A6, BO, B6. Top row. 
Two-dimensional histograms of radius and the adiabatic invariant K = T/n 2 ^ . Darker pixels represent a higher mass for a particular 
radius and K. This plot shows the wide variations in K at r > 100 pc, where cold flows and shock-heated gas coexist at a given radius. 
Radiative cooling allows the gas to cool and collapse in the center that accounts for the decrease in simulations A6 and B6. The material 
at r > 100 pc and K < 10 3 ' 3 K cm 2 corresponds to the cold flows inside filaments that illustrates that virialization occurs at different radii 
depending on its origin. Middle row. Two-dimensional slices of entropy. The circle denotes r200= 615 pc and 576 pc for simulation A and 
B, respectively. The virial shock exists at approximately r200 in the adiabatic models; however it shrinks to ~2/3 of r200 when we consider 
radiative cooling. Bottom row: Two-dimensional slices of number density of baryons. 



3.2. Variations in the virial shock 

Using the adiabatic invariant, K = T/n 2 ^ 3 , which we 
label "entropy" , allows us to differentiate between gas 
accreting from voids and filaments. As a precaution, we 
note that K is not an invariant when 7 varies; however, 
this is not the case in our simulations in which we permit 
molecular hydrogen cooling. Here molecular fractions re- 
main low, < 1CP 3 , and 7 w 5/3 even in the densest re- 
gions. The top row of Figure 3 depicts the variance of 
K with respect to radius in two-dimensional histograms, 
where the intensity of each pixel represents the mass hav- 
ing the corresponding K and r. The middle and bottom 
rows display two-dimensional slices of K and density, re- 
spectively, through the densest point in the halo. The 
virialized gas from the voids has low density and does 
not significantly contribute to the mass averaged radial 
profiles. Figure 3 illustrates this gas at r ~ r 2 oo and 
K > 10 4 ' 5 K cm 2 . The gas in filaments has lower entropy 
than the rest of the halo at r > 150 pc and K < 10 33 K 
cm 2 . In equation (2), the pressure in the surface term is 
the constant at a given radius. The accreting, denser, un- 
shocked gas in filaments has lower temperatures than the 
more diffuse accreting gas. The gas remains cool until it 
shocks and mixes well inside r v ; r and as small as ~7" v j r /4 
in the most massive filaments. Similar characteristics 



of cold accretion flows have been noted and discussed by 
Nagai & Kravtsov (2003a), Keres et al. (2005), and Dekel 
& Birnboim (2006). 

Entropy in the exterior of the halo differ little between 
adiabatic and cooling runs outside of r coo i. But as the gas 
falls within r coo i, it cools and condenses, which gives a 
lower entropy, and the r-K histograms and entropy slices 
display this clearly. Another significant difference in the 
cooling simulations is the contraction of the virial shock 
by a factor of 1 /3 when compared to adiabatic runs. This 
is caused by the contraction of the cooling gas. Here the 
cold filaments penetrate to even smaller radii. This is 
also evident in the radial density profiles of Figure 2. 

3.3. Virial Heating and Turbulence 

In order for a system to remain in virial equilibrium 
as it grows in mass, additional gravitational energy is 
balanced through two possible mechanisms: heating of 
the gas (£ > 0) and increasing the kinetic energy of 
the gas (T > 0). We differentiate between two main 
cases of virialization by comparing the cooling time, 
icooi = kT/nA, of the system to the heating time, iheat = 
T v i T /T v \x ~ |M vir /M vir in the case of rapid mass ac- 
cretion. Birnboim & Dekel (2003), Dekel & Birnboim 
(2006), and Wang & Abel (2007) find that radiative cool- 
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Fig. 4. — (a) simulation A. (b) simulation B. Radial (top) and tangential (bottom) velocity distributions of the most massive halo at z 
= z c jj2 (l e ft) and z = z c : t,ya (right). The heavy, blue lines are the distributions of the adiabatic models, and the light, black lines are 
from the radiative cooling models. Solid and dotted lines correspond to the velocity distributions of DM and baryons, respectively. These 
distributions can be decomposed into single or multiple Gaussians, depending on substructure. This demonstrates that violent relaxation 
occurs for the baryons as well as the DM. The narrower distributions of the baryons is due to the dissipation in shocks. 



ing rates are greater than heating rates from virialization 
for halos with masses below 10 12 M Q . 

1. Thermalization (t coo \ > theat) — When no efficient 
radiative cooling mechanisms (e.g. H2, Lya, He I) exist, 
the system virializes by injecting energy into heat £ and 
partly into kinetic energy T. In the process, the halo 
becomes pressure supported and virialized. Traditional 
galaxy formation scenarios only consider this thermaliza- 
tion while neglecting the kinetic energy term of equation 
(2). However it is important to regard kinetic energy, 
even in adiabatic models, as the gas violently relaxes. 
Turbulence velocities are similar to the velocity disper- 
sion of the system and contributes notably to the overall 
energy budget as seen in adiabatic cases in Figure 1. 

2. Turbulence generation (t coo \ < theat) — When a 
cooling mechanism becomes efficient, the system now 
dispenses its thermal energy and loses pressure support 
within r coo i. The gas will cool to a minimum equilib- 
rium temperature. As the cooling halo collapses and ra- 
dial velocities increase, the gas still lacks enough kinetic 
and thermal energy to match the gravitational energy 
and surface term in equation (2). The gas becomes more 
turbulent in order to virialize. We see this in the second 
row of Figure 1, where turbulent energies are significantly 
increased well inside the halo in the cooling models as 
compared to the adiabatic calculations. 

Through virial analyses, we have shown that turbu- 
lent energies are comparable, if not dominant, to ther- 
mal energies in galaxy formation. In the next Section, 
we further investigate the significance and nature of the 
turbulence through velocity distributions and decompo- 
sitions in order to study any small-scale anisotropics in 
the internal flows. 

4. VELOCITY DISTRIBUTIONS 

In CDM cosmogony, collisionless dark matter domi- 
nates the gravitational potential and oscillates as it sta- 
bilizes. Lynden-Bell (1967) showed how a collisionless 
system undergoes violent relaxation if embedded within a 
rapidly time- varying potential. Individual mass elements 
do not conserve energy during violent relaxation, only the 
entire system conserves energy. This behavior random- 
izes the energies of the mass elements, and statistical 



mechanics makes the resulting energy (velocity) distri- 
bution to tend to Maxwellian. Furthermore, the system 
"forgets" its original configuration during virialization or 
the incorporation of a lesser halo. Later studies have 
inferred two baryonic scenarios of virialization. First, 
violent relaxation and the accompanying phase mixing 
also applies to the gaseous component (van den Bosch et 
al. 2002; Sharma & Steinmetz 2005). In the other case, 
the gaseous component dissipates all turbulent motions 
and finally rigidly rotates as a solid body with a velocity 
appropriate to the overall spin parameter (e.g. Loeb & 
Rasio 1994; Mo et al. 1998; Bromm & Loeb 2003). 

Figure 4 shows the velocity distributions of the dark 
matter and baryonic components of the halo at z Cj H2 in 
the left column and z c ^ ya in the right column. It also 
overplots the simulations with adiabatic and radiative 
cooling. We plot the radial and tangential velocity dis- 
tribution on the top and bottom rows, respectively. The 
velocities are taken with respect to the bulk velocity of 
the halo. We also transform the velocity components to 
align the z-axis and total angular momentum vector of 
the DM halo. 

The radial velocity distributions at z — z Ci H2 are ap- 
proximately Maxwellian in both dark matter and gas 
with a skew toward infall. The infall distributions are 
shifted by ~1 km s^ 1 in the cooling case when com- 
pared to adiabatic. However at z = z c .Lya, the effects 
of Lya cooling become more prevalent in the halo when 
compared to H 2 cooling, shifting the radial velocity dis- 
tribution by ^5 km s -1 that is caused by faster infall. 
These distributions have two components that represent 
virialized gas and infalling gas in filaments. We further 
discuss this in the next Section. 

The tangential velocity distributions are nearly 
Maxwellian in all cases except for the dark matter in 
simulation A at z = z c ^ ya (right panels in Figure 4a). 
This deviation from Maxwellian arises from two major 
mergers that occur between 25 and 85 Myr (z = 17 
21) before the final collapse. The residual substructure 
from the major merger causes three distinct populations 
with Gaussian distributions centered at -0.2, +13.6, and 
-6.7 km s _1 with a = 11.6, 4.2 and 3.6 km s , respec- 
tively. These distributions clearly do not resemble a solid 
body rotator, whose velocity distribution would contain 
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Fig. 5. — Two-dimensional slices of velocity divergence (V • v) at z = z Cj Ly Q for simulations AO, A6, BO, and B6. The fields of view are 
1.49 and 1.69 kpc for simulations A and B, respectively. The circles mark the radius r2oo = 615 (576) pc for simulation A (B). Shocks are 
clearly denoted by large, negative convergence values. In the adiabatic cases, these shocks mainly exist at large radii where the gas from 
the voids and filaments virializes. When we consider radiative cooling, supersonic turbulence increases the frequency of shock fronts in the 
interior of the halo. 



all positive velocities. In other words the turbulent ve- 
locities exceed the typical rotational speeds. 

Distributions in dark matter are broader than the gas 
in both simulations and collapse redshifts as expected 
because for the gas we only give the bulk velocities and 
do not add the microscopic dispersion (cf. van den Bosch 
et al. 2002; Sharma & Steinmetz 2005). 

4.1. Halo and Filament Contrasts 

The dark matter velocity distributions are typical of 
a virialized system with the majority of the matter hav- 
ing a Maxwellian distribution with a dispersion corre- 
sponding to the main halo (Boylan-Kolchin & Ma 2004; 
Dieman et al. 2004; Kazantzidis et al. 2004). Substruc- 
ture appears as smaller, superposed Gaussians, which 
are stripped of its outer material as it orbits the parent 
halo. Dynamical friction acts on the substructure and 



decreases its pericenter over successive orbits, and the 
subhalo is gradually assimilated in the halo. 

The filaments penetrate deep into the halo and provide 
mostly radial infall inside r2oo ■ They do not experience a 
virial shock at r2oo, and this contrast is apparent in the 
radial velocity distributions. When we restrict our analy- 
sis scope to the filaments (i.e. r > 150 pc and K < 10 3 3 
K cm -2 ), the radial velocity distribution (Figure 4) is 
skewed toward infall, centered at -15 km/s, which is ap- 
proximately the circular velocity of the halo. The rest 
of the gas outside of this region in r — K space has al- 
ready been virialized, shock heated, and roughly exhibits 
a Maxwellian distribution, centered at zero, with its asso- 
ciated substructures. Hence the mass in filaments dom- 
inate the radial velocity distributions at negative values 
in Figure 4. 




Fig. 6. — Two-dimensional slices of the component perpendicular to the slice of vorticity (V X v) at z = z Cj L ya for simulations AO, A6, 
BO, and B6. The fields of view and circle radii are the same as in Figure 5. This quantity emphasizes the large- and small-scale turbulent 
eddies in the halo. 



4.2. Turbulence 

Radial inflows can create turbulence in the halo. Fil- 
aments provide an influx material with distinct angular 
momentum. This gas virializcs in the presence of an 
already turbulent medium that has a relatively high spe- 
cific angular momentum at r > r2oo/4. The Rayleigh 
inviscid instability criterion requires 

dj 2 

-rj- > for rotational stability, (9) 

where j is the specific angular momentum. If this is not 
satisfied, the system will become unstable to turbulence. 
The onset of turbulence can be delayed if viscosity were 
large enough so that Reynolds numbers are below the 
order of 10 2 or 10 3 . However there are many modes of 
instability if equation (9) is not met, and even a gas 
with low Reynolds number will eventually become fully 
turbulent (Shu 1992). 



Velocity dispersions can characterize the general mag- 
nitude of a turbulent medium, but its local nature is 
better detailed by applying the Cauchy-Stokes decompo- 
sition, 

u' = u+iwxh+i(V-u)h+ix>-h (10) 

that decomposes the velocity field into bulk motion u, 
vorticity w = V x u, expansion and contraction V • u, 
and a distortion V without change in volume. Here the 
vector h describes the separation between gas parcels at 
position x and x', and u' is the velocity at x'. 

We relate V • u, which is plotted in Figure 5, to a 
convergence timescale through the continuity equation, 

p + V(/ni) = 0, (11) 

that can be rewritten in terms of the total derivative 
D/Dt as 

I^ = -V.u. (12 > 
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Dp/Dt describes density changes along the fluid flow 
lines, and the 1/ p factor converts this change into an in- 
verse timescale on which local densities e-fold. We denote 
— (V • u) _1 as the "Lagrangian convergence timescale" 
(LCT). Converging flows (V • u < 0) arc ubiquitous 
within the halo. On large scales, the smallest LCTs on 
the order of 20 kyr exist at the virial shock, adjacent 
to both the filaments and voids. In the cooling models, 
these timescales are also small in turbulent shocks well 
within r 20 o- Analogous to the dynamical time, typical 
shocked LCTs decrease toward the center as density in- 
creases. 

The local magnitude and nature of turbulence is fur- 
ther illustrated by the vorticity lu, whose component per- 
pendicular to the slice, u> v , is shown in Figure 6. The lo- 
cal rotation period, 47r/|w| is also helpful to quantify and 
visualize the nature of the flow. The vorticity equation 
reads 

^+Vx(wxu) = — Vp x Vp, (13) 
at p A 

where the source term is non-zero when the density and 
pressure gradients are not aligned, i.e. baroclinic vortic- 
ity generation. This occurs at and near shocks through- 
out the halo, regardless of radiative cooling. In the adi- 
abatic models, vorticity exists even at modest but suf- 
ficient resolution in the large pressure-supported cores 
(see M.turb for r < 100 pc in Figure 1) and generates a 
turbulent medium with Mturb — 0.3. In the cooling mod- 
els, this large-scale vorticity is still present but increases 
in the collapsing core. As shocks become abundant in 
the center, we do not see any dampening of kinetic en- 
ergy. Perhaps this mechanism maintains turbulent mo- 
tions during virialization, even in the presence of dissi- 
pative shocks. In Figure 6, adjacent, antiparallcl fluid 
flows, i.e. a sign change in u y , are ubiquitous, which 
visually demonstrates that turbulence exists throughout 
the halo. The length scale of these eddies decrease with 
increasing density as with the LCTs. 

Hence we believe significant turbulence generated dur- 
ing virialization should be present in all cosmological ha- 
los. The cooling efficiency of the gas, the total halo mass, 
and partly the merger history determines the magnitude 
of turbulence. We discuss some implications of virial tur- 
bulence in the following section. 

5. DISCUSSION 

We have investigated the virialization of early pre- 
galactic cosmological halos in this paper with a suite 
of AMR simulations with varying chemistry and cooling 
models and collapse epochs. When analyzing the local 
virial equilibrium of the halo, we do not assume that it is 
in equilibrium but explicitly calculate all of the relevant 
terms in the virial theorem. In both adiabatic and ra- 
diative cooling cases, we find that the kinetic (turbulent) 
energy is comparable, if not dominant, with the ther- 
mal energy. Turbulence appreciates as radiative cooling 
becomes efficient because thermal energy alone cannot 
bring the system into virial equilibrium. In this case, the 
gas attempts to virialize by increasing and maintaining 
its kinetic energy. 

Besides violent relaxation, at least two other hydro- 
dynamic processes will augment virial turbulence. The 
first occurs when radial inflow interacts with the virial- 
ized gas. Due to the Rayleigh criterion, the high angular 



momentum gas creates an instability when it is deposited 
by filaments at small radii. The second happens when 
minor and major mergers create Kelvin-Hclmholtz in- 
stabilities and drives additional turbulence (e.g. Rickcr 
& Sarazin 2001; Takizawa 2005). Our results show that 
this turbulence is acting to achieve close to virial equi- 
librium at all stages during assembly and collapse. 

Virial turbulence may be most important in halos 
which can cool rapidly when compared to virial heating 
from mass accretion. Interestingly all halo masses below 
~10 12 M Q that can cool by Lya emission satisfy this con- 
dition (Birnboim & Dekel 2003; Dekel & Birnboim 2006; 
Wang & Abel 2007). 

Turbulence appears to mix angular momentum effi- 
ciently so that it redistributes to a radially increasing 
function, and thus only the lowest specific angular mo- 
mentum material sinks to the center. This segregation 
allows a collapse to proceed as it were self-similar and 
basically non-rotating. Similar results have been re- 
ported in cosmological simulations of collapses of the first 
stars (Abel et al. 2000, 2002: Yoshida et al. 2003, 2006b; 
O'Shca & Norman 2007). Here and in protogalactic col- 
lapses, the turbulent velocities become supersonic. One 
would expect even higher Mach numbers in larger poten- 
tial wells that still have i coo i < £d yn - 

The inclusion of the surface term allows us to study the 
virial equilibrium in the halo's interior where the grav- 
itational potential is not influential. Here our simula- 
tions show thermal and kinetic energies balancing the 
surface term and potential energy to achieve virial equi- 
librium. Before cooling is efficient, gas virially heats and 
its temperature can exceed the traditional virial temper- 
ature within r2oo/2 as seen in our adiabatic simulations. 
The consequences of this additional heating is substan- 
tial because halos can collapse and form stars before the 
virial temperature reaches the critical temperature, such 
as ^10,000K for Lya cooling. 

The timescale at which the center collapse occurs is 
crucial to the type of object that forms there. If the 
collapse occurs faster than a Kelvin-Hclmholtz time for 
a massive star (~ 300 kyr), a black hole might form 
from the lowest angular momentum gas. Conversely if 
the collapse is delayed by turbulent pressure, star forma- 
tion could occur in the density enhancements created by 
turbulent shocks. The ensuing radiative feedback may 
create outflows and thus slow further infall and possi- 
bly prevent the formation of a central black hole. The 
nature of the first galaxies poses an important question 
in the high-redshift structure formation, and to address 
this problem we must consider their progenitors - the 
first stars. 

5.1. Pop III Feedback 

Numerical simulations have shown that the first, 
metal-free (Pop III) stars form in isolation in its host 
halo. They are believed to have stellar masses ^100 
Mq (Abel et al. 2002; Omukai & Palla 2003; Tan & 
McKee 2004; Yoshida et al. 2006b) and produce -10 50 
photons s _1 that can ionize hydrogen and dissociate 
H 2 (Schaerer 2002). One-dimensional radiative hydro- 
dynamical calculations (Whalen et al. 2004; Kitayama 
et al. 2004) and recently three-dimensional radiative hy- 
drodynamical AMR (Abel et al. 2007) and SPH (Yoshida 
et al. 2006a) simulations found that pressure forces from 
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the radiatively heated gas drive a ~30 km s _1 shock out- 
wards and expels the majority of the gas in the host halo. 
Additionally the star ionizes the surrounding few kpc of 
the intergalactic medium (IGM). 

Pop III stellar feedback invalidates some of our as- 
sumptions in the calculations presented here, but the 
general aspect of kinetic energy being dominant should 
hold in the presence of these feedback processes. In a 
later paper, we will expand our simulations to include ra- 
diative feedback from primordial stars (cf. Yoshida et al. 
2006a; Abel et al. 2007) and the metal enrichment from 
pair instability supernovae (Barkat et al. 1967; Bond et 
al. 1984; Heger & Wooslcy 2002) of the IGM and subse- 
quent star formation. 

6. SUMMARY 

We have investigated the process of virialization in pre- 
galactic gas clouds in two cosmology AMR realizations. 
Our virial analyses included the kinetic (turbulent) en- 
ergies and surface pressures of the baryons in the sys- 
tem. The significance of each energy component of the 
gas varies with the effectiveness of the radiative cooling, 
which we quantify by performing each realization with 
adiabatic, hydrogen and helium, and H 2 cooling models. 
We highlight the following main results of this study as: 

1. Inside r coo i ; gas cannot virialize alone through heat- 
ing but must gain kinetic energy. It is up to a factor of 
five greater than thermal energy throughout the proto- 
galactic halos. This manifests itself in a faster bulk inflow 
and supersonic turbulent motions. 

2. In the radiative cooling models, supersonic turbu- 
lence (M = 1-3) leads to additional cooling within tur- 
bulent shocks. We expect turbulence in larger galaxies, 



up to 10 12 Mq, to be even more supersonic. 

3. Baryonic velocity distributions are Maxwellian that 
shows violent relaxation occurs for gas as well as dark 
matter. Turbulent velocities exceed typical rotational 
speeds, and these halos are only poorly modeled as solid 
body rotators. 

4. Virial shocks between the void-halo interface occur 
between r 20 o/2 and r 2 oo- Dense, cold flows in filaments 
do not shock-heat until well within r 2 oo and as small as 
r 2 oo/4. 

5. Turbulence generated during virialization mixes an- 
gular momentum so that it redistributes to a radially 
increasing function (the Rayleigh criterion). 

After the halo virializes, its central part will undergo 
turbulent collapse, such as in primordial star formation 
and galactic molecular clouds. These collapses should be 
ubiquitous in early structure formation as turbulence can 
be generated through virialization, merging, and angu- 
lar momentum segregation. We conclude that turbulence 
plays a key role in virialization and galaxy formation. 
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